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Abstract 

We consider the coagulation dynamics A + A — > A and the annihilation dynamics 
— * for particles moving subdiffusively in one dimension, both on a lattice and 
in a continuum. The analysis combines the "anomalous kinetics" and "anomalous 
diffusion" problems, each of which leads to interesting dynamics separately and to 
even more interesting dynamics in combination. We calculate both short-time and 
long-time concentrations, and compare and contrast the continuous and discrete 
cases. Our analysis is based on the fractional diffusion equation and its discrete 
analog. 



1 Introduction 

Diffusion-limited reactions in low-dimensional geometries have been studied 
intensely because they exhibit "anomalous kinetics," that is, behavior different 
from that predicted by the laws of mass action in well-stirred systems [1]. 
Among the simplest and most extensively studied are single species diffusion- 
limited coagulation (A + A — > A or A + A ^ A) [2] and annihilation (A + A — > 
0) [2,3]. These reactions, which show anomalous behavior in one dimension, 
are of particular theoretical interest because they lend themselves to exact 
solution in one dimension. Exact one-dimensional solutions can be obtained 
for diffusion in continuous [2,4-9] as well as in discrete [2,11-13] systems. A 
particularly elegant solution of the coagulation problem is provided by the 
method of intervals. This method focuses on the diffusive evolution of empty 
intervals, that is, of intervals that contain no particles [2,6-8,10]. The empty 
interval equation turns out to be linear and hence exactly solvable. Originally 
developed for continuous systems, the method has been extended to discrete 
lattices [11-13]. The method of intervals can not be directly adapted to the 
annihilation reaction, but a new formalism, the method of odd/even intervals 
which keeps track of the parity (even or odd) of the number of particles in 
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an interval, has recently been developed [8,9]. Again, the odd/even interval 
equation turns out to be linear and hence exactly solvable. 



The anomalies of the A + A problem in one dimension are typically displayed 
in two ways: one is through the time dependence of the reactant concentration 
c(t), which for the A + A — > A and the A + A — > reactions in infinite systems 
decays asymptotically as t^ 1 ^ 2 instead of the law of mass action decay The 
other is through the interparticle distribution function p(x,t), which is the 
(conditional) probability density for finding the nearest particle at a distance 
x on one side of a given particle. This function scales as x/t 1 ^ 2 , in typical 
diffusive fashion. In one dimension a gap develops around each particle that 
leads to a more ordered spatial distribution than the exponential distribution 
implicit in well-stirred systems and "explains" the relative slowing down of 
the reaction. 

In a parallel development, the problem of "anomalous diffusion" has also at- 
tracted a great deal of attention [14-16]. The universally accepted charac- 
terization of anomalous (as in "not ordinary") diffusion is through the mean 
squared displacement of a process x(t) for large t, 

Ordinary diffusion (a — 1, Ki = D) follows Gaussian statistics and Fick's 
second law for the probability density for finding the process at x at time t, 

lp<*t)-DlLp<?,t), (2) 

leading to linear growth of (x 2 (t)) with time. Anomalous diffusion is charac- 
terized by a nonlinear dependence. If < a < 1 the process is subdiffusive or 
dispersive; if a > 1 it is super diffusive. Anomalous diffusion is associated with 
many physical systems and is not due to any single universal cause, but it is 
certainly ubiquitous. Nor is anomalous diffusion modeled in a universal way; 
among the more successful approaches to the subdiffusive problem, which is 
the case we consider in this work, have been continuous time random walks 
with non-Poissonian waiting time distributions [14], and fractional dynamics 
approaches in which the diffusion equation (2) is replaced by the generalized 
diffusion equation [15-21] 

-P(x,t)=K a0 Dl-*—P(x,t) (3) 
where D}~ a is the Riemann-Liouville operator, 

o D}- a P(x, *) = _L|- /' dr -^^r , (4) 
* v ' ; T{a)dtJo (t-r) 1 ""' v ' 
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and K a is the generalized diffusion coefficient that appears in Eq. (f). Some 
limitations of this description and some connections between the generalized 
diffusion equation and continuous random walk formulations have been dis- 
cussed recently [22,23]. 

In this paper we consider a combination of these two phenomena, namely, the 
kinetics of A + A reactions of particles that move sub diffusively in one dimen- 
sion [24]. Some aspects of this problem have also been considered using the 
waiting time distribution approach [25,26]. Those solutions require approx- 
imations relating the reactant concentration to the distinct number of sites 
visited by a particle [25], or the waiting time distributions for single parti- 
cles to the waiting time distributions for relative motion [26]. Here we adapt 
the fractional dynamics approach to the problem and take advantage of the 
fact that the resulting generalized diffusion equations can be solved in closed 
form [24]. We consider both coagulation and annihilation reactions. 

In Section 2 we develop the fractional diffusion equations from which one 
obtains the concentration of reactant as a function of time. For the coagulation 
reaction we generalize the empty interval method [2,6] and, for the annihilation 
reaction, the odd/even interval method [8,9]. The solutions are presented in 
Section 3. Sections 2 and 3 deal with kinetics on continuum one-dimensional 
systems; in Section 4 we extend the analysis to reactions on a lattice. The 
differences between the continuum and discrete results are most pronounced 
at early times (long-time differences are due to finite size effects that get 
pushed to longer times if one deals with larger systems). We conclude with a 
recapitulation in Section 5. 



2 Evolution equations for subdiffusive particles 

2. 1 Empty interval method for coagulation reaction 

Consider first the coagulation reaction A + A — > A when the particles move 
by ordinary diffusion. The probability distribution function for the position x 
of any one A particle in the absence of reaction obeys the diffusion equation 
(2). The coagulation problem can be formulated in terms of the probability 
E(x, t) that an interval of length x is empty of particles at time t. This "empty 
interval" function obeys the diffusion equation [2,6] 

^E(x,t) = 2D^E(x,t). (5) 

The derivation of this equation is straightforward and recognizes that an empty 
interval is shortened or lengthened by movement of particles in and out at 
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either end according to the dynamics described by Eq. (2). 



It is instructive to recall this derivation as given in Refs. [6,7]. The reasoning is 
more transparent if the diffusion processes that lead to the change dE(x, t)/dt 
are described on a lattice with lattice spacing Ax and sites labeled by integers. 
In this lattice E n (t) gives the probability that sites I through n, for example, 
are empty at time t. Since the event E n+ \ contains the event E n , the proba- 
bility that sites 1 through n are empty but that there is a particle at site n+1 
is E n+ i — E n . The particles move randomly to the nearest lattice site with a 
hopping rate 2D /(Ax) 2 , the rate D /(Ax) 2 to the right being equal to the rate 
to the left. If the particle at n + 1 moves to site n in a short time interval At, 
E n will decrease by [D/(Ax) 2 ](E n — E n+ i)At. On the other hand, E n may 
increase by the departure of a particle at site n to site n + 1; the associated 
increase is [D / (Ax) 2 ](E n ^ 1 — E n )At. The same entry and exit processes can 
take place at the other end of the interval. Combining these four processes 
then gives for the change in E n 

The coagulation reaction fixes the boundary condition. Coagulation involves 
adjacent occupied sites and either of the two particles hopping onto the other, 
thereby disappearing. Two adjacent sites can both be occupied, or one can be 
occupied and the other empty (probability 2(Ei — E 2 ) because this can occur 
in two ways) or both can be empty (probability E 2 ). Since this covers all 
possibilities, the probability that two adjacent sites are occupied is therefore 
1 — 2Ei + E 2 . The change in E 1 is thus 

and for this equation to fit the pattern of Eq. (6) it is necessary to require the 
boundary condition 

E = 0. (8) 

Letting x = nAx and both Ax and At — > then leads to the continuous 
equation (5) with the boundary condition 

E(0,t) = 0. (9) 

The other boundary condition requires that the population of particles be 
nonvanishing: 

E(oo,t) = 0. (10) 

The empty interval dynamics is thus essentially the same as that of individual 
particles in the absence of reaction, but with different boundary conditions 
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and a diffusion coefficient 2D that reffects the fact that the relative motion of 
two diffusive particles involves twice the diffusion coefficient of each particle 
alone. 



The connection between the empty interval function and the observables of 
interest is obtained as follows [6,7]. Back in the discrete formulation, the prob- 
ability that a site is occupied (i. e. not empty) is 1 — E\, and the concentration 
of particles is therefore 



Ax ' dx 



(11) 

x=0 



where we have exhibited the continuum limit. The interparticle distribution 
function p n (t) is the probability that the nearest neighbor to one side of a 
given particle is n lattice spacings away. This probability is related to E n by 



oo oo 



E n = cAx T,Pm- (12) 

k=n+l m=k 

The first sum insures that the next neighbor is at least k sites away and the 
second that k is at least n+ 1. The normalization takes into account that the 
average distance between the particles is the reciprocal of the concentration, 

OO J 

(nAx) = ^2 np n Ax = -. (13) 

n=l C 

Equation (12) can be inverted: 

(E n+1 -2E n + E n _ l ) d 2 E(x,t) 

CPn= Ax "~ ' < t )P( X > t )= dx 2 ( M ) 

where again we have exhibited the continuum limit. 

The construction of the kinetic equation for E(x,t) for subdiffusive particles 
proceeds along arguments analogous to those used in the diffusive case. Again, 
one follows the motion of the particles in and out at the ends of the empty 
interval according to the dynamics (3). However, for particles that move sub- 
diffusively, the concept of a hopping rate (steps per unit time) is not defined 
because the number of steps performed up to time t by a subdiffusive particle 
goes as t a with a < 1 [23]. Therefore, the derivation of the empty interval 
equation must be adjusted accordingly. Proceeding as above in the discrete 
formulation, we see that the probability that sites 1 through n are empty but 
that there is a particle at site n + 1 is still E n+ i — E n . However, the rate 
of decrease of E n due to a particle that moves from site n + 1 to site n is 
now described by the generalized diffusion process implicit in Eq. (3), that is, 
[K a Dl~ a / (Ax) 2 ](E n — E n+ i)At. Collecting arrival and departure processes 
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at both ends of the interval then leads to the equation 

(E n+1 -2E n + E n _ 1 ). (15) 



AE n 2K a Dj- a 



At (Ax) 2 
In the continuum limit this then becomes 

^E(x, t) = 2K a o Dl~ a -^E(x, t). (16) 

The boundary conditions are exactly as before, Eqs. (9) and (10), as are the 
relations (11) and (14) to the observables. 



2.2 Odd/ 'even interval method for annihilation reaction 



The empty interval method can not be applied to the annihilation reaction 
A + A — > because annihilation leads to a discontinuous growth of empty 
intervals. However, recently a new method of odd/even intervals has been 
introduced that leads to exact solution in the diffusion-limited case. It is based 
on the construction of an equation for r(x, t), the probability that an arbitrary 
interval of length x contains an even number of particles at time t [8,9]. This 
probability does not change if two particles inside the interval react since this 
process does not change the even/odd parity of the interval. Because r(x,t) 
changes only by the movement of particles in or out of the ends of the interval, 
arguments similar to those that lead to Eq. (5) lead to the same equation 
for r(x,t). This procedure can again be directly extended to the subdiffusive 
problem, and r(x,t) satisfies the same fractional diffusion equation as E(x,t): 

^r(x,t) = 2K a0 Dl-<»-^r(x,t). (17) 

One boundary condition is the same as for E(x,t) and is obtained via the 
same arguments used earlier: 



r(0,t) = l. (18) 

The other boundary condition in general differs from that of the empty interval 
method and depends on the initial condition. In particular, it is determined 
by the parity of total initial number of particles. This parity never changes. A 
random initial placement leads to an equal probability that the system forever 
contain an even or and odd number of particles, so that 

r(oo,f) = ^. (19) 
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The concentration of particles is related to r(x,t) precisely as in Eq. (11): 

dr(x, t) 



c(t) = 



dx 



(20) 



x=0 



3 Solution of the fractional diffusion equations and subdiffusive 
reaction kinetics 



3.1 Coagulation reaction 



The solution of Eq. (5) with boundary conditions (9) and (10) can be written 
as [18] 

roc rt(x / \f2K a )~ 2 / a 

E(x, t) = J dy [W(x -y,t)- W(-x - y, t)} E(y, 0) + J ' " uj a/2 { V )d V , 



where 



W(x,t) 
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\x\ 


(1-f.f) 


V2K a t<* 


(0,1) 



(21) 



(22) 



up(x) 



px 
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1//5) 



(23) 



and H is the Fox H-function [15,16,18,27]. Taking into account that in Laplace 
transform space (indicated by a tilde over the function and the letter u as 
variable) 



W(x,u) 



u 



a/2-1 



cxp 



u 



a/2 



X 



(24) 



and co/3(u) = exp(— u 13 ), the solution E(x,t) in Laplace space adopts a much 
simpler form: 

____ o rCO , \ ] 

E(x, u) = — dy ( e -\ x -y\ s - e-\ x+y \ s ) E(y, 0) + - exp (-xs) (25) 
2u Jo v ' u 



where 



a/2 



S = 



V2Ka" 



(26) 
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From Eqs. (11) and (25) one finds 



c(u) = - l-sE(s, 0) 
u 1 



(27) 



where E(v,t) = / °° dxe vx E(x,t) is the the spatial Laplace transform of 
E(x,t). Equivalently, from Eq. (14) one obtains 



c(u) = £[l-p( a ,0)] 



(28) 



where A = c(0) and p(u, 0) is the spatial Laplace transform of the initial 
interparticle distribution function p(x, 0). 

A commonly considered initial interparticle distribution is the random (Pois- 
son) distribution of average concentration A, for which E(x, 0) = e~~ Xx and 
Po(x) = \e~ Xx . For this initial distribution either Eq. (27) or Eq. (28) can be 
used to obtain 



c(u) 



A 



(29) 



and c(t) is given in closed form in terms of the Mittag-Leffler function [16,28] 
of parameter a/2: 



c(t) = \E a/2 (-X V^W 2 ) 



(30) 



The Mittag-Leffler function can be calculated by the series expansion [16,28] 



j^ o r(fca + l) 



(31) 



When a — 1 one recovers the usual result for diffusion-limited coagulation [4,5] 
since the Mittag-Leffler function of parameter 1/2 is E 1 / 2 (—x) = exp(x 2 )erfc(a;). 

Another initial particle distribution of interest is an ordered or periodic ar- 
rangement, p(x, 0) = 5(x — 1/A). In this case Eq. (28) leads to 



c[u) 



A 

u 



1 — exp 



u 



a/2 



and c(t) is given in terms of the Fox H-function H]_'® [18,27,7]: 



C(0 . rrlQ 

— = 1-^11 



(32) 



1 


(l.f)" 


W2K a t a 


(0,1) _ 



(33) 
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For normal (a = 1) diffusive particles one recovers the usual result for diffusion- 



limited coagulation [5] since H 





(M)" 


z 






(0,1) _ 



erfc(f). 



It is useful to exhibit explicitly the long-time and short-time behavior of the 
particle concentration. The long-time behavior can be extracted via Taube- 
rian theorems from the small--u behavior independent of initial condition and 
therefore in particular common to both of the initial conditions considered 
above, 

? ,-<*/2-i 

*"> ~ 7B?r (34) 

from which it follows that 

--a/2 

c(t) ~ ; ? r-. (35) 

v^r(i-f) 

For diffusive particles (a = 1) one recovers the well-known result c(t) ~ t^ 1 / 2 . 
The short-time behavior of the particle concentration is, not surprisingly, 
markedly different for the two initial conditions. For the random initial con- 
dition an expansion of Eq. (30) using (31) gives 

Notice that this implies a infinite initial reaction rate because dc(t)/dt — > oo 
for t — >■ 0, a reflection of the fact that a random distribution of particles in- 
cludes a large probability of proximate particles. On the other hand, using the 
asymptotic expression of the Fox H-function H\'1 for large argument [18,27] 
in the solution for the periodic initial distribution of particles, Eq. (33), one 
finds 



a I a 

.2 



a/(2-a) 

z 2/(2-a) 



+ . . . 

(37) 



where z = [\(2K a t a ) l / 2 Y l . In this case dc(t)/dt — > for t — > 0, since initially 
there are no proximate particles. 

The curves in Fig. 1 shows the time evolution of the concentration c(t) for 
the coagulation reaction and the two initial conditions considered here. The 
differences at early times are evident, as is the coagulation of the curves into 
a single asymptotic form at long times. 
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o 



i < 1 1 1 < r 




Fig. 1. Logarithm of the survival probability c(t)/X for coagulation dynamics of 
subdiffusive particles with a = 1/2, K a = 1, and initial concentration A = 1/20, 
versus the logarithm of time. The initial arrangements of the particles are Poissonian 
(solid line and circles) and periodic (broken line and squares). The lines are the 
continuum solutions given by Eqs. (30) and (33). The symbols are the on-lattice 
solutions for a finite lattice with L = 200 and Ax = 1 with periodic boundary 
conditions. 

The time-asymptotic form of the empty interval function can be used to obtain 
the asymptotic behavior of the interparticle distribution function via Eq. (14). 
From Eq. (25) one deduces that 



d^ EM = Yu 



I" dy ( e - lx - yls - e - lx+yls ) E(y, 0) + —e- x ', (38) 
Jo v ' u 



which for u — > reduces to 



d 2 ~ s 2 



(39) 



From Eq. (14) and the fact that the inverse transform of (39) can be expressed 
in terms of a Fox H-function, one finds 



1 



c(t) p(x,t) ~ -H™ 



X 


(Lf)" 


V2K a t«/ 2 


(2,1)_ 



(40) 



which in combination with Eq. (35) immediately leads to 



X 


(l-«,f) 


V2K a t«/ 2 


(0,1) 



(41) 



10 




1.0 



0.0 



0.0 



0.5 



1.0 



1.5 



2.0 



2.5 



3.0 



z 



Fig. 2. Long-time scaled interparticle distribution function for the coagulation re- 
action for several values of the anomalous diffusion exponent. Proceeding upward 
from lowest to highest curves along the y axis intersection: a = 1, 0.95, 0.8, 0.5, 0.2. 
Note that the distribution for a = 0.2 on this scale is nearly indistinguishable from 
the completely random distribution exp(— z) (dotted curve). 

for t — > oo. This result is independent of the initial distribution of particles. 
We can also write this in the scaled form 



where z = c(t)x is the scaled interparticle distance and p a (z)dz = p(x,t)dx. 
This stationary form is shown in Fig. 2 for several values of the diffusion 
exponent a. 

The interparticle distribution function conveys the interesting "anomalies" of 
the problem most clearly. For a random distribution of particles on a line this 
distribution is exponential. In particular, the most probable interparticle gaps 
are the smallest. For diffusion-limited reactions on a line it is well known that 
the scaled distribution deviates in two ways from the exponential behavior. 
First, a gap develops around each particle, and the distribution vanishes near 
the origin (see a = 1 curve in the figure), indicating an "effective repulsion" of 
particles. Second, the probability of large gaps decays much more rapidly than 
exponentially: the decay goes as a power of exp(— z 2 /2). In the subdiffusive 
case decreasing a leads to the diminution of the gap around each particle, 
that is, to a weakening of the effective repulsion and to a behavior that ap- 
pears closer to that of a random distribution in the short-interparticle-distance 
behavior. This is evident in the progression of the curves with decreasing a 
shown in the figure. Furthermore, the probability of large gaps decays as a 




(42) 
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power of exp(— x 2 ^ 2 ~ a ^), thus neatly interpolating between the purely ran- 
dom exponential decay state as a — > (since then p(x,t) — > c(t)e~°^ x ) and 
the more ordered state corresponding to diffusive particles at a — 1. 



3.1? Annihilation reaction 



It is straightforward to deduce that the solution r(x, t) of Eq. (17) with initial 
condition r(x,0) — \ + E(x,0)/2 and boundary conditions (18) and (19) is 
simply r(x, t) = | + E(x, t)/2 where E(x, t) is the solution found in the previ- 
ous section. We can therefore here rewrite Eq. (27) for the Laplace transform 
of the concentration of particle undergoing the annihilation reaction as 

c(u) = -[l-sr(sM (43) 
u 

where r(v,t) is the the spatial Laplace transform of r(x,t). 

For a random initial distribution of mean concentration A the initial distribu- 
tion is r(x, 0) = \ + \e~ 2Xx (note that this goes to 1/2 as x — > oo since there 
is an equal probability that we start with en even or odd number of particles, 
and it goes to as x — > 1 since zero particles is an even number). Hence 

~° {U) = M + 2A^W2 » (44) 
which can be inverted analytically to yield 

c(f) = XE a/2 (-2\j2K~ a t a l 2 ^ . (45) 

Comparing this result with Eq. (30) for the coagulation reaction we see that 
the rate of change of the concentration for annihilation is exactly twice that 
for coagulation at all times, a result that has been noted for ordinary diffu- 
sion [8,9]. 

For a periodic initial distribution with separation 1/A between particles we 
have 



r(x, 0) 



2j + l- \x 2j/X <x<(2j + l)/\ 

(46) 

{-2j-l + Xx (2j + l)/X<x<(2j + 2)/X 



This is a periodic function of period 2/A, so that [29] 



*») = ^^ (47) 
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Fig. 3. Logarithm of the survival probability c(t)/X for annihilation dynamics of 
subdiffusive particles with a = 1/2, K a = 1, and initial concentration A = 1/20, 
versus the logarithm of time. The initial arrangements of the particles are Poisso- 
nian (solid line and circles) and periodic (broken line and squares). The lines are the 
continuum solutions given by Eq. (45) and by the numerical inversion of Eq. (48). 
The symbols are the on-lattice solutions for a finite lattice with L = 220 (filled sym- 
bols) and L = 200 (open symbols) with periodic boundary conditions and Ax = 1. 
Note that with these values of L and A the initial number of particles and hence 
the total number of particles is odd (filled symbols) or even (open symbols) at all 
times. 

Carrying out this integral and substituting into Eq. (43) immediately leads to 



where we have written the second equality to make comparison with Eq. (32) 
more apparent. We have not been able to invert this expression analytically 
to find c(t) in closed form. 

Figure 3 shows the temporal evolution of c(t) for the two initial conditions 
when a = 1/2. The evolution for a periodic initial distribution has been ob- 
tained by numerically inverting Eq. (48). 

The long-time behavior is again seen to be independent of initial condition and 
can be obtained directly from Eq. (43). Since f(s, 0) — > (2s)~ x when s — > 0, it 
follows that as u — > 0, 




(48) 




-a/2-1 



(49) 
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which in turn implies that 



t -a/2 

C(t) 7 r-, (50) 

2v^r(i-f)' 

when t — > oo. Comparison with Eq. (35) shows that at long times the con- 
centration for the annihilation reaction is exactly half that of the coagulation 
reaction. This is in agreement with the general result obtained above com- 
paring the rates of change of the concentrations. The short-time behaviors for 
the two initial conditions are different, as already evident in Fig. 3. For the 
random initial condition an expansion of Eq. (45) leads to 



[compare with Eq. (36)]. For the periodic initial condition expansion of the 
denominator of Eq. (48) for large u and comparison with Eq. (32) leads to 



cm 2 ^y/<°-y„- 2)exp 



a 



z 2/(2-a) 



+ 

(52) 



where again z = [XfiKj ) 1 / 2 ]- 1 [compare with Eq. (37)]. 



4 Reactions in a lattice 



It is sometimes convenient not to carry out the continuous limit in the formula- 
tion of the reaction kinetics for coagulation or annihilation reactions. Not only 
may the actual physical system be discrete, but simulations usually involve 
discrete lattices, and finite reactant size effects (i.e., small distance scale ef- 
fects) are more appropriately dealt with through discrete formulations [10-13]. 
Indeed, a discrete viewpoint was the starting point of the continuum equations 
considered in the previous section. Furthermore, the discrete formulation al- 
lows consideration of different reaction rules such as A$A — > AAA [12]. 

The hierarchy of differential- difference equations on a lattice of L sites is 



-^ = a D l t - a [I n+1 (t) - 2I n (t) + J n _i(f)] 1 < n < L - 1 (53) 

for n = 1,2,..., where a = 2K a /(Ax) 2 , I n (t) = E(x n = nAx,t) for the 
coagulation reaction, and I n (t) = r(x n = nAx, t) for the annihilation reaction. 
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The other boundary condition, at n — L, depends on which reaction is under 
consideration and will be stated below. The discrete form [6] 

was already introduced in Eq. (11). From here on, we take Ax = 1. 



4-1 Coagulation in a segment 



The problem of the coagulation dynamics in a discrete segment has been 
solved for the normal diffusive case [10-13]. Here we present the solution for the 
coagulation process of subdiffusive particles on a chain of L sites with periodic 
boundary conditions, which is described by the hierarchy (53) together with 
the additional boundary condition 

hit) = (55) 

We solve this set of equations by means of the ansatz 

In(t)=J2 b n(u)Ea(-aUJt a ) (56) 

U! 

which is closely related to the method of separation of variables for the subd- 
iffusive fractional differential equation [16,20]. Again, E a is the Mittag-Leffler 
function of parameter a. Taking into account that 

f t E a (-Pn = -(3 Dl- a E a (-pt a ), (57) 

the resulting eigenvalue problem can be solved in the same way as for normal 
diffusion [12,11]. Notice that E a (x) = e x for a = 1 so that the ansatz (56) 
takes the well-known form for normal diffusive particles (for which a = 1). 
The full solution is 

I n (t) = (l - I) + An sin (m|(n - L)) E a {-auu m t) (58) 

with 

u m = 2[1 - cos(m7r/L)]. (59) 
The coefficients A m are determined from the initial conditions. 
For an initially random distribution of particles of concentration A we have 
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7 n (0) = (1 - A) n for n = 0, 1, • • • , L - 1, so that 

2 (l-A)W7r 2 + (-irL 2 ln 2 (l-A) 
mlJ " mTT m 2 7r 2 + L 2 ln 2 (l-A) lbUj 

For an initially periodic distribution of particles separated by 1/A lattice sites 
we have J„(0) = 1 — \n for \n < 1 and J„(0) = for An > 1. In this case 

= +1 2ALBi, l( ,^/AL) 

771 7T 



Figure 1 shows the decay of the concentration c(i) given by Eq. (54) with (58) 
when the initial distribution of subdiffusive particles on the lattice is random 
and when it is periodic. At longer times both solutions approach those of the 
continuum equations and would do so more clearly if we were to exhibit the 
solutions to our difference equations for larger lattices. The marked deviation 
of the lattice solution from the continuous asymptotic t^ a ^ 2 decay occurs when 
the concentration on the lattice approaches its asymptotic (minimum) value of 
one particle on the entire lattice, c(t — > oo) = 1/L. This time can be estimated 
using the continuum solution to occur at time t x such that c(t x ) = 1/L. In the 
figure this separation of solutions clearly does occur when c/A = 1/XL — 0.1, 
i.e., when log 10 (c/A) = —1.0. 

While the long-time behavior of the discrete and continuous solutions coin- 
cide (independently of initial condition) until the reactant concentration is 
extremely low, the short-time behavior of the two solutions is expected to 
differ, especially when the initial distribution of particles is random. The dif- 
ference has been established for particles undergoing normal diffusion [13]. To 
determine the short-time behavior on the lattice we assume that the lattice 
is initially full, so that I n (0) = for n > 1. (The reasoning follows along the 
same lines for an arbitrary concentration of particles, but the equations are 
more cumbersome.) Taking the time Laplace transform of Eq. (53) gives 

I n+1 (u) - (2 + — )I n (u) + /„_!(«) = 0, n > 1. (62) 
a 

This set of difference equations must be solved with the boundary conditions 
Io(u) = l/u and Il(u) = 0. Since we are interested in the short-time ( large-it ) 
behavior, we can set L — > oo without appreciably affecting the outcome. With 
these boundary conditions we essentially follow the procedure of Ref. [13] to 
find 

In particular, ii(it) = a/u 1+a + . . . for u — > oo, so that with the help of 
appropriate Tauberian theorems we find that I\{t) = at a /T(l + a) + . . . for 



2 + 



it 



'u 2a u a \ 1/2 



(63) 
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t — > 0. Using Eq. (54), we calculate the short-time behavior of the reactant 
concentration 

Note that with an initially full lattice there is no distinction between an ini- 
tially random and initially periodic distribution of particles. For an arbitrary 
random initial concentration one still finds that c(t)/X — 1 ~ t a . Comparison 
with Eq. (36) shows that at short times, as in the case of normal diffusion [13], 
the concentration on a lattice decays more slowly than in the continuum. This 
is not easily visible on the scale of Fig. 1. 



4-2 Annihilation in a segment 



Now consider the annihilation reaction A + A — > and suppose that initially 
there are N particles on a segment of length L. The hierarchy (53) is now 
augmented with the second boundary condition 

!1, N = even 
(65) 
0, N = odd 

where we note that now I n (t) = r(x n = nAx,t). This condition is of course a 
result of the fact that the reaction does not change the parity of the number 
of surviving particles. We proceed as in the coagulation problem and write the 
solution as 



4(0 = In + A ™ sin rn-(n - L) E a (-au m t) (66) 

m=0 V L 1 



where 



1, N = even, 

I* n = I n (oo) = { (67) 
1-2 N = odd, 

oj m is given in Eq. (59), and the A m are coefficients determined from the initial 
conditions. 

For a random initial condition of particle concentration A we have I n (0) = 
\ + |(1 - 2A) n for n = 0, 1, • • • , L - 1, so that 



A m (t) = 



(-1) N - (1 - 2A) L 1 m 2 7r 2 + [l - (-l) m+N ] L 2 ln 2 (l - 2A) 



mn m 2 7r 2 + L 2 ln 2 (l -2A) 

(68) 
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For a periodic initial distribution with particle separation 1/ A we have 

'2,7 + 1 -An 2j/X < n < (2j + 1)/A 
-2j - 1 + An (2j + 1)/A < n < (2j + 2)/A 



4(0) = 



(69) 



which leads to 



A m (t) = (-if 1 - (-i) N (-ir 



2 XL ( rrm \ 
ri^ tan \2X~L 



(70) 



Figure 3 shows the evolution of the concentration as given by Eq. (54) with 
(66) for a random initial distribution (circles) and a periodic initial distribution 
(squares). A number of points are noteworthy. First, the solution of course 
decays to a finite constant (one particle remaining) when there is initially an 
odd number of particles in the system, and to zero when the initial number of 
particles is even. The deviation of the solutions on a finite segment from those 
in an infinite continuum again set in around the time when c = 1/L (strictly 
speaking around 2/L for the even case). Note that in the continuum solution 
one does not distinguish between the even and odd particle number cases: the 
boundary condition (19) is an average of the two, and so the continuum curves 
should be compared to the average of the discrete ones for the two boundary 
conditions. 



5 Conclusions 



We have considered the coagulation dynamics A + A^A and the annihilation 
dynamics A+A — > for particles moving sub diffusively in one dimension. This 
scenario combines the "anomalous kinetics" and "anomalous diffusion" prob- 
lems, each of which leads to interesting dynamics separately and to even more 
interesting dynamics in combination. The fractional diffusion equation plays 
a central role in our analysis and allows the exact calculation of the density 
c(t) within this formulation. We have calculated c(t) explicitly for all times for 
the coagulation reaction in a one-dimensional continuum for two initial dis- 
tributions, a random (Poisson) distribution and a periodic distribution. This 
calculation is based on the empty interval method. Using the odd/even inter- 
val formulation, we have also obtained an explicit solution for the annihilation 
reaction with a random initial distribution. For a periodic initial distribution 
we are only able to calculate the Laplace transform of c(t) analytically and 
must perform the inversion numerically. 

For the coagulation reaction the empty interval method also leads to an ex- 
plicit solution for the interparticle distribution function p(x,t). At long times 
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we find a universal expression for this distribution (i.e., one independent of 
initial distribution), in terms of a single scaled variable. Anomalous diffusion 
is characterized by the exponent a introduced in Eq. (1), ordinary diffusion 
corresponding to a = 1. Deviations from ordinary diffusion lead to a curi- 
ous interplay. On the one hand, with decreasing a (and hence increasingly 
subdiffusive motion) the decay of the particle density towards extinction be- 
comes increasingly slower and in this sense increasingly different from law of 
mass action behavior. On the other hand, for the coagulation reaction the 
spatial distribution of initially randomly distributed reactants remains more 
Poissonian for all time as a decreases; indeed, as a deviates from unity the 
relatively empty regions around each particle that are tantamount to (and 
indeed explain) anomalous kinetics in the usual diffusion-limited case become 
more populated. 

We have also considered the problem of subdiffusion-limited reaction on dis- 
crete lattices. We have presented a hierarchy of differential-difference equations 
for the empty interval probability I n (t) = E(x n ,t) and the even/odd interval 
probability I n {t) = r(x n ,t), and have expressed the solutions I n (t) as a super- 
position of subdiffusive modes whose decay is governed by the Mittag-Lemer 
function in the same way that the exponential function governs the decay of 
ordinary diffusive modes. As in ordinary diffusion [13], the associated concen- 
trations c(t) differ from the continuum solutions at early times. At long times 
the lattice solutions deviate from those of the continuum systems when there 
is only a small number (one or two) of particles remaining in the system. 
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